
% script3 plots the cluster for different angles for one chosen frequency
% the angles start from +80 and goes until -80
% (this code is very slow!!)

clear all; close all;clc;
        rad = 0.1;
cc=0;
thre=0.01;
load hrir_final021.mat
freqs = [200 300 500 600 800 1000 1200 1500 2000 3000];
nturns=[0.1 0.4 0.8 0.7 0.8 1.2 1.3 1.4 1.8 2.2];
% % freqs= [ 400 600 800 1000 1200]; %just a random frequency
% freqs= [ 100]; %just a random frequency
% nturns=[0.5 0.9 2.1 4.05]
% nturns=[0.65 0.9 1.2 1.4 1.5]
% nturns=[1.42]
% titles = {'A. Female Speech','B. Male Speech','C. Music'};

fs=44.1e3;
% Nump='ABCDEFG'
soundwave=randn(fs*15,1); % BB

f1=figure('position',[100 100 800 800]);
f2=figure;

G=[-80 -65 -55 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 55 65 80];
Elevation_number = 9;
for ifreql = 1:length(freqs)
    freq = freqs(ifreql)
    tdur = 40/freq;
    for i=1:25
        i
        H2l=squeeze(hrir_l(i,Elevation_number,:));  H2r=squeeze(hrir_r(i,Elevation_number,:));
        x1=filter(H2l,1,soundwave);
        x2=filter(H2r,1,soundwave);
        % % % %      short time fourier transform:
        [SL,SR,freq_pos, flag01]=mySTFT(x1,x2,fs,thre,tdur);
        %      normalization:
        [SL1,SR1]=Normalize(SL,SR);
        [Y Ind] = min(abs(freq_pos - freq));
        X_P = [real(SR1(Ind,flag01>0))' imag(SR1(Ind,flag01>0))'];

        figure(f1)
        subplot(2,5,ifreql)
        plot(X_P(:,1),X_P(:,2),'.','Markersize',6);
        [IDX, C] = kmeans(X_P, 1,'Distance','cityblock','Replicates',5);
        F_centers(i,:) = C;
        set(gca,'Xtick',-1:0.5:1)
        set(gca,'Ytick',-1:0.5:1)
        % title([Nump(ifreql) '. ' num2str(freq), ' Hz'],'fontsize',8);
        title([num2str(freq), ' Hz'],'fontsize',8);
        %     axis([-1 1 -1 1])
        hold on;
        plot(-1:1,[0,0,0],'k', 'LineWidth', 0.5);
        plot([0,0,0],-1:1,'k','LineWidth',0.5)
        axis square
        box off
        txy = C;
        flag_plot = 1;
        if ifreql<3 & mod(i,2)==0
            flag_plot = 0;
        end
        if flag_plot==1
            s=text(txy(1)-0.02,txy(2)-0.12,[num2str(G(i)) '^{\circ}'],'horizontalalignment','center');
        end
        %     s=text(txy(1),txy(2),[num2str(G(i)) '^{\circ}'],'horizontalalignment','center')
        set(s,'fontsize',10)
        grid on
        set(gca,'Xtick',-1:0.5:1)
        set(gca,'Ytick',-1:0.5:1,'fontsize',10)
        axis([-1.1 1.1 -1.1 1.1])
        xlabel('Real','fontsize',13); ylabel('Imaginary','fontsize',13);

        pause(0.01)
    end

      
    nturn=nturns(ifreql);
    %==============================
    Norder=6;
     out=f_get_spiral(F_centers,Norder,nturn);
     x_spiral=out.x_spiral;
     y_spiral=out.y_spiral;
    figure(f1)
    subplot(2,5,ifreql)
    plot(x_spiral,y_spiral,'color',[0.2 0.2 0.2]) ; % nturns crossings, including end point
%     stop0
    hold on
    exportgraphics(gcf, 'CC_Figure_5_Various_Freq.png', 'Resolution', 300,'ContentType', 'vector')

%     figure(f2)
%     subplot(2,2,ifreql)
%     bar(r_nonlinear)
%     hold on;plot(xx_fine,fitted_r,'r')
    pause(0.1)
    clear x1 x2

    % eval(['save Freq_centers_hardware' num2str(freq) ' F_centers x_spiral y_spiral'])
end
%     if i_soundtype==1
%         exportgraphics(gcf, 'Female_model.png', 'Resolution', 300,'ContentType', 'vector')
%     elseif i_soundtype==2
%         exportgraphics(gcf, 'Male_model.png', 'Resolution', 300,'ContentType', 'vector')
%     else
%         exportgraphics(gcf, 'Music_model.png', 'Resolution', 300,'ContentType', 'vector')
%     end
